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ABSTRACT 

We present a novel method for the light-curve characterization of Pan-STARRSl Medium Deep 
Survey (PS1 MDS) extragalactic sources into stochastic variables (SV) and burst-like (BL) transients, 
using multi-band image-differencing time-series data. We select detections in difference images associ¬ 
ated with galaxy hosts using a star/galaxy catalog extracted from the deep PS1 MDS stacked images, 
and adopt a maximum a posteriori formulation to model their difference-flux time-series in four Pan- 
STARRSl photometric bands gpi,rpi,ipi, and zp\. We use three deterministic light-curve models 
to fit burst-like transients; a Gaussian, a Gamma distribution, and an analytic supernova (SN) model, 
and one stochastic light curve model, the Ornstein-Uhlenbeck process, in order to fit variability that 
is characteristic of active galactic nuclei (AGN). We assess the quality of fit of the models band-wise 
source-wise, using their estimated leave-out-one cross-validation likelihoods and corrected Akaike in¬ 
formation criteria. We then apply a K-means clustering algorithm on these statistics, to determine 
the source classification in each band. The final source classification is derived as a combination of the 
individual filter classifications, resulting in two measures of classification quality, from the averages 
across the photometric filters of 1) the classifications determined from the closest K-means cluster 
centers, and 2) the square distances from the clustering centers in the K-means clustering spaces. For 
a verification set of active galactic nuclei and supernovae, we show that SV and BL occupy distinct 
regions in the plane constituted by these measures. We use our clustering method to characterize 
4361 extragalactic image difference detected sources in the first 2.5 years of the PS1 MDS, into 1529 
BL, and 2262 SV, with a purity of 95.00% for AGN, and 90.97% for SN based on our verification 
sets. We combine our light-curve classifications with their nuclear or off-nuclear host galaxy offsets, to 
define a robust photometric sample of 1233 active galactic nuclei and 812 supernovae. With these two 
samples, we characterize their variability and host galaxy properties, and identify simple photometric 
priors that would enable their real-time identification in future wide-held synoptic surveys. 

Subject headings: methods: statistical, techniques: photometric, (galaxies:) quasars: general 


1. INTRODUCTION 

With the advent of the LSST era, human-intervened 
classification (with the exception of citizen science) will 
become untenable, requiring automated source identifi¬ 
cation in large volumes of data in archival catalogs, as 
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well as in real-time data. Recent increases in computa¬ 
tional resource availability and efficiency have enabled 
near-com plete automated transient discovery in large 
surveys (Bloom et al. 2012). Machine-learning methods 
are slowly replacing human judgement for transient clas- 
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that use a number of photometric and non-photometric 
features have been demonstrated to achieve classifica¬ 
tions with very high purity. 

As the number of detected transients grows very large 
in wide-held time domain surveys, complete spectro¬ 
scopic follow up becomes impossible due to limited re- 
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sources and faint magnitude limits. Classification meth¬ 
ods using time-series data alone are favorable, and have 
been app lied in_the past to a broad range of sources; 


Choi et al. 


(2013) discuss the identification of AGN 
via damped-random wal k parameteriza t ion of imagc- 
differencing light curves, Andrae et al. (2013) on the 
applicability of single and multiple (Jr nstein-Uhlenbeck 


(OU) processes to AGN light curves, and Butler & Bloom 
(2011) on the separation of AGN from variable stars in 
photometric surveys through damped-ra ndom walk pa¬ 


ramet erization. For supernovae (SNe), [Kessler et al. 
(2010) discusses various photometric template fitting 
methods that enable their identification with particular 
SN classes. 


Bailer-Jones (2012) gives a powerful method to pick 
suitable models for deterministic or stochastic lightcurves 
using leave-one-out cross-validation. This method can 
be used to determine source type, based on the likeli¬ 
hoods of the various models for a given lightcurve. So 
far, the applicability of selecting sources based on time- 
series modeling h as bee n limited to single-band detec¬ 
tions (Choi et al. 2013 ), or have typically u sed magni¬ 
tude time-series data (Butler & Bloom 20111, which are 
undefined for negative difference-fluxes. Also, computa¬ 
tional limitations typically have lead to the use of only 
single models as predictors for class, or only using simple 
statistical criteria for model assessment, rendering clas¬ 
sification schemes prone to the possibility of systematic 
misclassification. 

We attempt here to use time-series methods alone to 
classify sources into general categories of burst-like (BL) 
or stochastically variable (SV). These light-curve classes 
capture the variability behavior of the two most common 
extragalactic sources detected in image-differencing sur¬ 
veys, AGN and SNe. We present a novel method that 
separates BL and SV sources with high purity, using su¬ 
pervised machine-learning methods. Using multi-band 
difference-flux in the gpi,rpi,*pi and zp± bands, we se¬ 
lect BL and SV from 4361 difference-image sources with 
galaxy hosts. For lightcurves in each band, we estimate 
the fitnesses of analytical BL models relative to the OU 
process, using both their estimated leave-out-one cross- 
validation likelihoods (LOOCV), and corrected-Akaike 
information criteria (AICc). We show that the use of sim¬ 
ple analytical models with suitably chosen priors, which 
mimic the approximate shapes of BL light curves (pre¬ 
dominantly SNe), is sufficient for segregating them from 
SV, thereby obviating the need for exact models depen¬ 
dent on specific BL subclasses. The model statistical 
characterizations are combined across sources using a 


K-means clustering algorithm (Kanungo et al. 


2002), to 


provide robust source classifications in each filter. The 
filter-wise classifications are then averaged to give final 
source classifications. Based on our BL and SV pho¬ 
tometric classifications and host galaxy offset cuts, we 
define a photometric sample of SNe and AGN, which we 
then use to define observational priors for future surveys 
such as LSST. 

This paper is organized as follows: In f[2] we discuss 
the details of the survey and our data pipeline; in m we 
discuss time-series models used to describe BL and SV 
difference-flux light curves; in )j4] we elaborate on com¬ 
puting LOOCV and AICc to estimate model fitness; in 
we discuss our classification method, characterize the 


TABLE 1 

Pan-STARRSI Medium Deep 
Survey Field Centers 



J2000 

J2000 

Field 

RA 

Declination 


Degrees 

Degrees 


MD01 

02 h 23 m 30 s 

-04 deg 15' 

MD02 

03 h 32 m 24 s 

-27 deg 48' 

MD03 

08 h 42 m 22 s 

44 deg 19' 

MD04 

10 h 00 m 00 s 

02 deg 12' 

MD05 

10 h 47 m 40 s 

58 deg 04' 

MD06 

12 h 20 m 30 s 

47 deg 07' 

MD07 

14 h 14 m 48 s 

53 deg 04' 

MD08 

16 h ll m 08 s 

54 deg 57' 

MD09 

22 h 16 m 45 s 

00 deg 16' 

MD10 

23 h 29 m 14 s 

00 deg 25' 


properties of extragalactic variables and transients, and 
combine our light-curve classifications with host galaxy 
offsets in order to define a robust photometric sample 
of SNe and AGN; and in (J6]we describe the source and 
host galaxy properties of our variability/offset selected 
SN and AGN samples, and report on priors that can aid 
in their real-time classification in future surveys. 

2. THE PAN-STARRSI SURVEY AND TRANSIENT ALERT 
DATABASE 



55000 55200 55400 55600 55800 56000 
MJD 


Fig. 1.— The Pan-STARRSI survey has a staggered 3-day 
cadence in the gpi,rpi,ipi, and zp\ bands corresponding to 6 
observations per month per filter, while yp\ is observed during 
bright-time. The observations we use in this paper extend from 
2009 September 14 till 2011 November 17. In this paper yp\ is not 
used due to the relatively sparse cadence as compared to the other 
filters. 


Th e Pan-STARRSI (PS1) telescope (Hodapp et al. 


20041 is a 1.8 meter diameter telescope on the summit 


of Haleakala, Hawaii with a //4.4 primary mirror, and 
a 0.9 m secondary, delivering an image with a diameter 
of 3.3 degrees onto 60, 4800 x 4800 pixel detec tors, with 


10/Ltm pixels that subten d 0.258” each (Tonry & Onaka 


2009||Hodapp et al.|2004|). The observations are obtained 
itt 


through a set of 5 broadband filters gpi,r’pi,*pi, 2 pi,ypi, 
each with a limiting magnitude per nightly epoch of 
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~ 23.5 mag. Although the filter system for PS1 has much 
in common with that used in previous surveys, such as 
the SDSS, there are substantial differences. For more 


tech ni cal de tails refer to Stubbs et al. (2010), and Tonry 
eTaLl(p0T2). 


The PS'1 survey has two operating modes, 1) the 37r 
survey which covers 37 t square degrees at <5 > —30 de¬ 
grees in 5 bands with a cadence of 2 observations per 
filter in a 6 month period, and 2) the Medium Deep Sur¬ 
vey (MDS) which obtains deeper multi-epoch images in 
5 bands of 10 fields, each 8 square degrees, listed in Ta¬ 
ble [lj designed for both extensive temporal coverage, and 


once the alert is registered, forced photometry can de¬ 
tect both positive and negative fluxes. These criteria 
remove the majority of “bogus” detections due to non- 
astrophysical sources, such as camera defects, cosmic 
rays, and image-differencing artifacts. The PS1 alerts 
are published to an online alerts database located in 
Harvard. Our automated pipeline then downloads the 
alerts database to our local database servers at Uni¬ 
versity of Maryland on a nightly basis. The alerts are 
then processed and additional value added measurements 
are made on the data to enable easy characterization 
of sources via a SQL-IDL-C++ pipeline (Fig. [2]). The 


weather, the accessible fields are observed with a stag- 

multi-band deep-stack catalogs ( 

Heinis et al. 

in prep.|) 

gered 3-day cadence in each band during dark and gray 

to derive host associations. Other statistics sue 

i as color 


time (gpi,rpi on the first day, ipi on the second day, 
zpi on the third day, and then repeat with <?pi, rpi), and 
in the y^i band during bright time. On average, the ca¬ 
dence (Fig. [T]) is 6 observations per filter per month, with 
a 1 week gap during bright time, during which time the 
Medium Deep fields are observed exclusively in ypi. We 
require the dense time-series (cadence«few days) for ro¬ 
bust variability-based classifications, thereby making the 
MD survey our survey of choice. 

The PS1 MD data are processed using the image pro¬ 
cessing pipeline (IPP) located in Hawaii. The IPP per¬ 
forms flat-fielding and detrending on each of the indi¬ 
vidual images using white light flat-field images from a 
dome screen, in combination with an illumination cor¬ 
rection obtained by rastering sources across the field of 
view. Bad pixel masks are applied, and carried forward 
for use in the stacking st age. A fter det ermin ing an ini¬ 


tial astrometric solution (Magnier et al. 2008), the flat- 
fielded images are then warped onto the tangent plane 
of the sky, using a flux conserving algorithm. The im¬ 
age scale of the warped images is 0.250 arcsec/pixel. In 
the MD fields, all images from a given night are col¬ 
lected with eight dithers. This allows the removal of de¬ 
fects like cosmic rays or satellite streaks, before they are 
combined into a nightly stack using a variance-weighted 
scheme. Nightly stacks of images, each with a 8 square 
degree field of view, as well as seasonal deep stack ref¬ 
erence images are created, which are then transferred 
to the Harvard Faculty of Arts and Sciences Odyssey 
Research Computing cluster, where they are processed 
through a frame subtraction analysis using the photpipe 
image differencing pipeline originally d evel o ped fo r the 


evolution, and higher moments of magnitude and flux 
are also computed and stored in our database. Web 
pages that derive custom cuts on the data based on host 
properties, host offsets, color, magnitude, and time vari¬ 
ability properties are also updated nightly. Our custom 
query page can be used to query the database and display 
column-wise sortable results on a web page. The page 
can also be used to visualize the data in our database 
using simple 2 dimensional plots or histograms that are 
created in IDL which are displayed on a web page. Fi¬ 
nally, the transient alerts are classified based on their 
light curves using our time-series method discussed in 


2.1. Pre-Processing the Alerts for Classification 

The two most commonly occurring extragalactic time- 
varying sources are SNe and AGN, which fall under the 
broad categories of burst-like and stochastically vary¬ 
ing, respectively. However, these broad classifications, in 
combination with host galaxy offsets, also enable us to 
discover more rare and exotic variables and transients. 
For example, nucl ear BL sources should include tidal 
disruption events (Gezari et al. 2012 Chornock et al. 
20141, off-nuclea r BL sources may contain gamma-ray 
burst afterglows (|Cenko et al.||2013| |Singer et al.||2013|) 
and off-nuclear SV sources may be offset AGJN from a 
post- merger recoiling supermassive black holes (Blecha 
etaL|2011 |. 

We identify extragalactic alerts by cross-matching the 
18, 058 alerts detected in the first 2.5 years of the PS1 
MDS with galaxies detected in our custom multi-band 


deep-stack star/galaxy catalogs (|Heinis et al. 


m prep 


Garg et al. 12007 

). Significant flux excursions are then de- 

built from CFHT’s u-band and the five PSl bands. The 

tected in tne cli 

Terence images using 

a modified version 

detection threshold, defined by the yfi distribution, is 

of the DoPhot ( 

Schechter et al. 

1993 

Rest et al. 

2013) 

equivalent to a SNR of 1.9<r. The photometry is then 


photometry package, and they are tagged as a source, if performed using SExtractor (|Bertin & Arnoutsj|1996|) in 
they satisfy the following conditions: Kron elliptical apertures which are also used m cross 

matching alerts with the objects in the catalog. The cat¬ 
alog contains ss 10 7 objects which have been classified as 
stars or galaxies with over 90% accuracy for sources with 
magnitudes < 24 mag, using an optimized SVM classifi¬ 
cation scheme that takes into acc ount the shape, co lor 
and magnitudes of the detections ( Heinis et al.|in prep. 


• Positive detections with a signal-to-noise ratio 
(SNR) > 5 in at least three images within a time 
window of 15 days. 

• Detections in at least two filters. 

• No previous alert at that position. 

While for source detection, DoPhot requires at least 3 
consecutive positive detections within a 15 day period, 


Thus we only select alerts with ipi -host < 24 mag, where 
the star/galaxy classification is reliable. We identified 
4361 extragalactic alerts using the catalog, which we then 
characterized as SV or BL using multi-band difference- 
flux time-series. Note that we do not include extragalac- 
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Fig. 2.— The PS1-UMD data pipeline. The data are relayed from the IPP via the photpipe pipeline that provides transient alerts, as well 
as performs forced photometry and image differencing. At UMD the data are downloaded, enhanced with statistical parameterizations, and 
assimilated into SQL databases using a C++ framework, which are then queried using IDL or PHP for interactive web-based analysis. 


tic alerts with unresolved hosts in our sample, such as 
quasars, or “hostless” alerts, those with either a faint 
host galaxy (*pi_h os t > 24 mag), or located outside the 
elliptical region that defines their host galaxy. 

To characterize the sources, we model their difference- 
flux time-series obtained from forced photometry in each 
of the gpi,rpi,ipi, and zpi bands, and then combine the 
characterizations across the filters. Forced photometry, 
done by the photpipe pipeline, is obtained by performing 
PSF fitting photometry on all difference images in each 
band, at the location of any transient candidate, in order 
to fully exploit all available difference-flux time-series in¬ 
formation on the alert, using data from prior to the alert 
detection. In our method, we decided to use difference- 
fluxes instead of differential magnitudes because stochas¬ 
tic light curves can have negative difference-fluxes (if they 
were brighter in the reference image), for which AB mag¬ 
nitudes cannot be defined. 

Before we perform the classification on the difference- 
flux light curves, we pre-process them to remove arti¬ 
facts, as well as to make them conform with SV and BL 
model priors. Many difference-flux light curves also con¬ 
tain singular large difference flux outliers which affect 
model classification significantly. To remove them, each 
difference-flux point yi in a given light curve is compared 
to its previous and next difference-flux values, pi-\ and 
Ui+ 1 , and their photometric errors dpi-\ and dj/j+i,with 
the criterion that at least one of \y t — Pi~\\ < I0dyi-i, 
or | yi — Pi+i\ < 10dpi+i is satisfied for the difference- 
flux point to be accepted as non erroneous. Since most 
difference-flux errors are much larger than this cutoff and 
most differences in successive difference-flux values are 
much smaller than this, we can ensure that the light 
curve is unaffected, while the outliers are removed. Since 


the starting and ending points of light curves cannot be 
subject to one of these criteria, we discard them after 
removing the erroneous difference flux points. Also, our 
condition will not weed out the turn-on of BL lightcurves, 
because only \pi — yi-i\ < 10dyi-i may not be satisfied 
due to a short rise-time, however, |pi — yi+i\ < lQdpi+i 
will be satisfied for all BL lightcurves, given the Pan- 
STARRS1 medium-deep survey observations repeat once 
every 3 days, in the gp±, rpi, ipi, zpi bands. Finally, we 
also transform the light curves such that the minimum 
difference-flux is 0. This is done so as to make them 
conform to the limits for the priors for BL light curves, 
especially for light curves where the BL source was active 
in the reference image. 

In our analysis we only use light curves which have at 
least n = 20 distinct difference-flux measurements post¬ 
processing in any filter, so as to ensure that this is at least 
four times as large as the maximum number of parame¬ 
ters k m ax used in any of the models (Maximum number 
of parameters in any time-series model (jj3]) is 5). This 
is done to prevent over-fitting of the data, which may 
result in model comparisons not being meaningful. Also, 
n = 20 is not a restrictive limit for classification pur¬ 
poses since this is a factor ss 2 smaller than the average 
number of photometric measurements in any filter for 
all 4361 sources, which is ss 36. Fig. [3 shows the his¬ 
togram of number of distinct points in the ffpi, rpi, zpi, 
or zpi filters for all the PS1 transient alerts associated 
with galaxy hosts that pass our cuts. In the next section 
we discuss the time-series models that we use to classify 
the difference-flux light curves. 

3. TIME-SERIES MODELS 

Since our goal is to classify extragalactic time-varying 
sources into two broad classes, BL or SV, we assess the 
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Fig. 3. — Sources plotted in the figure satisfy our criteria for classification: a light curve with n > 20 points in all 4 filters. MD10 has 
no points which satisfy our criteria due to only a single full season of coverage in the first two and a half years of the PS1 MDS. 


general shapes of the light curves by comparing their 
similarities to SN-like bursting behavior, or to AGN-like 
damped-random-walk type behavior. While fitting an 
exact model involves a large number of parameters which 
may be unknown, and may necessitate a large number of 
data points, the general shape of a BL light curve can 
be approximated to certain simpler analytical functional 
forms (Gaussian, Gamma distribution, and generic ana¬ 
lytic SN model); and that of an SV lig ht curve ap prox¬ 


imated by an OU p rocess (Uhlenbeck & Omstein 1930 
Andrae et al. 2013) as described in Table |T| In these 
models, we have ignored the effects of cosmological red- 
shift corrections and dust extinction. However this is 
acceptable since our goal here is to use the models only 
to distinguish between coherent single-burst type behav¬ 
ior from stochastic variability, while not assuming any 
underlying physical processes for the sources. Since the 
model parameters span several orders of magnitude, their 
priors are chosen to be uniform in the logarithm. 

The Gaussian is the simplest model that attempts to 
model the overall flux from a BL source as the sum of 
a constant background a, and bursting behavior charac¬ 
terized by a Gaussian with amplitude /3, center /i, and 
width cr. This however, does not account for the asym¬ 
metry in SN light curves; for example Type la SNe are 
better approximated by a sharp rise f r j se ~ 15 days (Gai- 


followed by a relatively slow decline in the flux in any 
band (t f a u ~ 30 days). To resolve this, we employ a 
Gamma distribution which is robust in modeling such 
light curves (Fig. |4|, reflecting varying degrees of asym¬ 
metry depending on the shape k and scale D parameters 
of the distribution. Another model is the Analytic-SN 
model, that uses distinct exponential rise and decline 
timescales, t v j Bfi and tf a n, and is particula rly well suited 


to mo deling non-la type SN light curves (Kessler et al. 


2010), although it is generic in its application. Despite 


the non-specific nature of these models, we find that their 
simple statistical descriptions of the difference-flux light 
curves of BL sources are sufficient in distinguishing them 
from AGN with low contamination. The use of three dis¬ 
tinct analytical BL models allows for a broader range of 
BL light-curve shapes, and is comparable to using inde¬ 
pendent statistical descriptions of the light curve through 
distinct parameterizations. Also, since the BL models 
are compared with the SV model, only their relative fit¬ 
nesses in describing the data are important. Should the 
necessity arise of classifying the objects into particular 
sub-classes of the broader SV-BL distinction, or that of 
extracting particular details about the paramet ers of a 
source light curve, exact mode ls for the sources (Kessler] 


et al. 2010 |Kelly et al. 2010) must be included in the 
comparisons, which although it is beyond the scope of 


better approxima ted by a snarp rise t r j se ~ to day s (Uai- comparisons, wmen aitnougn it is beyo 

tan|2013 jGaneshalingam et al.|20lT Hayden et al. 201 (If this paper, is a direct natural extension. 
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Fig. 4.— A range of BL light curve shapes can be modeled using 
Gamma distributions by varying the shape and scale parameters 
k, D. This is particularly applicable to the asymmetric rise and fall 
time-series patterns of SN light curves. 

The fluctuating behavior of optical li ght curves of AGN 
is well described by an OU process (Kelly et al .|2009 ), a 
first-order continuous-time auto-regressive process. The 
process can be described in terms of a driving noise field, 
parameterized by c, the square of the amplitude of the 


OU noise field, and a damping timescale r (Bailer-Jones 


2012). Mathematically, the evolution of the state vari¬ 
able Z(t) of the OU process is given by the differential 
equation 


dZ{t) = c 1/2 dW{t) - -(Z{t)-b)dt (1) 

T 

where W is a Wiener process, and b is the mean value of 
the process. For simul ating the OU pr oces s itself, we use 
the prescriptions from Bailer-Jones (2012). The method 
uses posterior analysis to improve the estimate of the 
state variable continuously, by using the observed flux 
y k - 1 at time-step t k ~ 1 to compute the posterior distri¬ 
bution of the state variable z k - i, which is subsequently 
used to update Zk- The OU process being a Gaussian- 
Markov process, Z(t) is characterized by Gaussian prob¬ 
ability distribution function G(/j,(Z(t)),V(Z(t))) where 
/x, V are the mean and standard deviation of the state 
variable at time t. 

We also determine whether the light curves are well 
fitted by a constant model that is representative of white 
noise. In the event that none of the light curve models 
is significantly better than white noise, the light curves 
are assumed to not pertain to any of the stochastic or 
bursting categories, and are classified as No-model (NM) 
sources. Figs. Hi and [7] show examples of SN, AGN, 
and NM classified sources and all the model fits. In each 
case, the best models are chosen based on robust sta¬ 
tistical criteria, and the final source class decided using 
a clustering machine learning scheme described in the 
following sections. 


4. MODEL LIKELIHOOD AND FITNESS ESTIMATION 

For all time-series models, including for the OU pro¬ 
cess, we assume a Gaussian error model to compute the 


Fig. 5.— An SN difference-flux light curve is reasonably well 
fit by all the models, but the BL models have higher LOOCV and 
lower AICc as compared to the OU process resulting in the light 
curve being classified as BL. 



Fig. 6.— Example of an AGN light curve that is well fit by the 
OU process and poorly by the BL models. 


model likelihoods. Although for an OU process, th e ac- 
tual likelihood is computed differently from this (|Bailer^ 
Jones||2012|), using only the photometric errors to com¬ 
pute the likelihood for all models, is justified, since the 
intent is to determine the model-mean that best mimics 
the light curve shape. The probability P{ykWk,d n ) of 
observing a difference flux yk, assuming Gaussian errors, 
is given by 


logP(y k \<T k ,0 n ) = log (—\ ^ Vk) ^ 

\oW27r/ 2a k 

where f k , yk, and a k are the model difference-flux, the 
observed difference-flux, and the standard deviation es¬ 
timates of the kt h data point. For the OU process we 
use n(Z(tk )) or the mean light curve, in the place of f k , 
to evaluate its likelihood. 


























7 


TABLE 1 

Difference-flux models 


Model 

Type 


Equation 

Parameters 

Prior Distributions 


Gaussian 

BL 

Flux(t) = a + /3e f 4 F)~ !° 2 

a 

p 

p 

(7 

Ulog( 0,10 5 )(flux counts) 
Ulog( 10 3 ,10 8 )(flux counts) 
Ulog(- 10 2 ,10 4 )(days) 
Ulog( 1,10 4 )(days) 

Gamma distribution 

BL 

Flux(t ) = 


a 

p 

p 

D 

k 

Ulog(0, 10 5 )(flux counts) 
Ulog(0 , 10 !, )(flux counts) 
Ulog(—10 2 , 10 4 )(days) 
Ulog(l, 10 2 ) 

Ulog(l, 10 2 ) 

a + P D k r(k ) 

Analytic-SN model 

BL 

Flux(t) 

-“ + ^l + e-(*-*o)/V„ e 

a 

P 

to 

t fall 
trise 

Ulog(0 ,10 5 )(flux counts) 
Ulog(0 ,10 9 )(flux counts) 
Ulog(—10 2 , 10 4 )(days) 
Ulog(l, 10 3 )(days) 
Ulog(l , 10 3 )(days) 

OU process 

SV 

dZ(t) = —^Z(t)dt + c 1 / 2 iV(t; 0, dt) 
where Z here is flux count. 

T 

C 

b 

MZ) 

V(Z) 

Ulog( 1,10 6 )(days) 

Ulog( 0,10 14 )(flux counts 2 ) 
Ulog( 0,10 8 )(flux counts) 
Ulog(0 ,10 8 )(flux counts) 
Ulog( 0,10 14 )(flux counts 2 ) 

No-Model 

White Noise 


Flux(t) = C 

c 

Ulog{ 0,10 8 )(flux counts) 



Time (days) 


Fig. 7.— Example of a difference-flux light curve that has a large 
number of image-differencing errors resulting in it being best fit by 
the No-Model. 


To assess the fitness of the models, we esti mate t heir 


corrected Akaike information criteria (AICc) (Burnham 

T 


& Anderson 12002 ) and leave-out-one cross-validation 


likelihoods (LOOCV) (Bailer-Jones 2012) over the 
difference-flux data for each source, filter-wise. The AIC 


(Eq, 


mode 


3) is a quantification of the information lost when a 


is used to represent a dataset. The AIC penalizes 
the maximum model log-likelihood InC by a factor that 
depends on the number of model parameters fc, thereby 
accounting for over-parameterization of the dataset. 


AIC = 2k- 2 InC 


(3) 


AICc = AIC + 


2 k(k + 1) 
n — k — 1 


(4) 


The AICc is a correction to the AIC, that corrects for 
the finite size of the dataset n relative to the number 
of model parameters k. Note, that models that better 
represent the dataset have smaller AICc values. The 
LOOCV, another independent measure of model fitness, 
is a measure of how well each difference-flux value can 
be predicted using the remaining difference-flux data and 
hence, is a more robust statistical measure of model fit¬ 
ness as compared to the AICc. Re-stated, the cross- 
validation can also be said to measure the predictive abil¬ 
ity of a mo del on a g i ven datas et, and is a nearly unbiased 
estimator (Kohavi et al.|1995) of the mean squared error 
for new observations. The LOOCV, more specifically, is 
the product of the piece-wise probability of obtaining in¬ 
dividual difference-flux measurements using a time-series 
model, while sampling the parameters from the posterior 
constituted by the model over the remaining points in 
the time-series. In LOOCV estimation of a model over a 
dataset yk containing K points, the likelihood L of the 
fcth data point is given by 


^ n—N 

Lk = P(y k \y-k,cr,r]) « — ^2 p (ykWk,0 n ,il) ( 5 ) 

n—1 

where rj is the time-series model, at is the error es¬ 
timate at each point, and 9 n are the model parameters 
drawn in the nth iteration from the Markov chain Monte- 
Carlo sampling of the posterior probability distribution 
of the model over the other K — 1 data points denoted 
by y_fc. The LOOCV of the model can then be obtained 
by multiplying the partition probabilities 



























k—K 


LOOCV= n Lk 


( 6 ) 


k —1 


Since the AICc and the LOOCV are measured inde¬ 
pendently of each other, they can be used simultaneously 
to assess model likelihood, thereby reinforcing model fit¬ 
ness assessment. The LOOCV for each model is esti¬ 
mated using a Markov chain Monte-Carlo (MCMC) using 
a standard Metropolis- Hastings a lgorithm to sample the 
posterior distributions ( jHastings 1970). The model pa¬ 
rameters are sampled from known distributions and the 
posterior probability LiPi is evaluated, where pi is the 
prior probability and Li is the model likelihood in the 
ith iteration. Parameters for the i + 1th iteration are ac¬ 
cepted with probability (Li+ipi+i) / (LiPi), failing which 
the parameters from the ith iteration are retained. We 
use a log-normal sampling distribution with a diagonal 
covariance matrix, with a ^ = 10~ 4 uniformly across all 
parameters, and all models. We find that this choice of 
a constant variance of the sampling distribution leads to 
stable cross-validation likelihood values. In accordance 
with log-normal sampling requirements, the parameter 
distributions defined in Tableware transformed between 
—oo,oo using a sigmoidal transform. 

The prior distributions for the BL and SV model pa¬ 
rameters can be assumed to be uniform in the logarithm, 
as we have in our simulations, or can be obtained by sam¬ 
pling the parameters at the posterior maxima for the BL 
and SV models, for known SNe (BL) and AGN (SV) 
training sets. The latter is advantageous if the entire 
set of sources is well represented by the training set, in 
that the number of iterations to convergence would be 
significantly reduced. However, we did not make this as¬ 
sumption in order to allow for the classification of BL 
and SV light curve types that may not occur in the veri¬ 
fication set, and only took care to ensure that the limits 
on the parameter ranges subsumed the parameter values 
that could occur in the dataset. 

Since the initial guesses for the model parameters in 
the MCMC may be far from the actual solution, a burn- 
in of 1000 iterations is employed for all model assess¬ 
ments. We determined that a large number of burn-in 
iterations is important to ensure sampling near the peaks 
of the posterior distribution, and is particularly impor¬ 
tant while using prior parameter distributions that are 
uniform in the logarithm, as we have done here. We 
determined that 10000 post burn-in iterations were suf¬ 
ficient for good model-fit convergence, after replicating 
the results with 2000 burn-in iterations, and 20000 post 
burn-in iterations. To ensure that the Markov chains 
were well mixed, we further comp uted the Gelman-Rubin 
statistic (Gelman & Rubin 1992) over 10 parallel chains, 
for each parameter in each model, and ensured that their 
values were less than 1.1. The calculation of the LOOCV 
is tedious and computationally expensive, and required 
us to parallelize our codes over a 300 core multi-node 
cluster. In addition our codes were written ground-up in 
C-|—|- and optimized for quick run-time. The classifica¬ 
tion of ss 7000 sources with ss 40 difference-flux points 
in each of the four bands, required ~ 4 hours. 


5. CLASSIFICATION METHOD 


Once the fitnesses of the models are estimated filter- 
wise on the difference-fluxes using the AICc and the 
LOOCV, we obtain two parameters per model for 5 time- 
series models in each of the four filters gp i, 7 'pi, ipi, and 
zp\. First, we remove the NM best fit sources by com¬ 
paring their model statistics with those of the BL and 
SV models. To do this we construct a relative sign vec¬ 
tor RVij for each source, in each filter, using the AICc 
and the logarithms of the LOOCVs (which we designate 
as LLOOCV): 


RV,f = { 

Sgn(LLOOCVGa U ssian,iJ LLOOCVnM,iJ ); 
Sgn(LLOOCVGamma,i,f - LLOOCV N M,ij), 
Sgn{LLOOC\^Analytic,i,f )) 

Sgn(AICc G aussian,iJ AIC Cn M i 

Sgn(AICcGamma,iJ AICcp[M,iJ ) 5 
sgn( AICc-Analytic,ij AICCj\[ M ,i,f 
sgn(LLOOCVou,ij ~ LLOOCVnm^j), 
sgn(AICcou,ij — AICcnm,ij) 

} (7) 


where i is the object id, / is the filter, and sgn 
denotes the sign function, defined to be +1 for posi¬ 
tive values and —1 for negative values. Ideally, for a 
BL source, RVbl = {+1, +1, +1, —1, —1, —1, ±1, 4=1} 
since the BL models will have a larger LLOOCV, and 
a smaller AICc when compared to the same for the NM, 
while for an SV source the relative sign vector should 
be RV sv = {±l,±l,±l,=Fl,Tl,=Fl,+l,-l}, i.e., the 
OU process better describes the light curve as com¬ 
pared to any of the BL models or the SV models. For 
sources where the NM is the best model, RVnm = 
{— 1) — lj~ Ij+Ij+Ij+I) - 1) + !}■ 

We then compute RVij for all sources, which are ag¬ 
gregated in filter-wise and fed into a K-means clustering 
supervised-machine-learning algorithm, using the num¬ 
ber of centers K = 3 in a swap method th at is repeated 
over 100 iterations (Kanungo et al. 2002 1 . The cluster¬ 
ing algorithm partitions the sources in the 8 -dimensional 
RVij space, into Voronoi cells to determine the centers 
of the distributions for BL, SV, NM, by minimizing the 
sum of squares of the distances of points Xi within cluster 
S m from the means of the clusters p m that correspond 
to the different classes of sources: 


k 

H \\ X l-Vm\\ 2 (8) 

l—l Xi eSm 

Each source is then assigned a class Cij as BL, SV, or 
NM depending on the center it is clustered around. The 
squared-distance of each source point i in filter /, Dij = 
\xi—p,cj\ 2 from the clustering center pej is a measure of 
how reliably it is classified as the particular type C , with 
a distance of Dij = 0 being the best, and larger distances 
indicating less reliable classifications. C\j and Dij are 
computed for each source, in each of the gpi,rpi,ipi, 
and zp\ bands independently. We choose to classify the 
sources filter-wise, and not using the statistical measures 
from all the filters at once, for the following reasons. 











9 


1. The behavior of each type of source, across the fil¬ 
ters, cannot be assumed to be uniform and hence, 
the clustering centers may differ significantly, 

2. Clustering in some filters may be more noisy than 
others, resulting in most sources being classified 
as no-model sources, thereby making these bands 
less favorable for classification purposes. For these 
filters, the no-model center would be repeated in 
place of an SV or a BL center. 

3. Some filters may be less noisy and show clustering 
only around two centers corresponding to BL and 
SV. Combining these filters with the noisier ones 
results in both, more uncertainty in clustering clas¬ 
sification (larger Di) and a larger number of mis- 
classifications. This is because, the uncertainty in 
the clustering classification caused by one or more 
filters confounds the otherwise clear classifications 
from the others. As a result, the clustering cen¬ 
ters are poorly determined in the joint parameter 
space of statistical parameters from all the filters. 
By performing the clustering in each filter sepa¬ 
rately, the classifications can be reinforced if they 
show agreement across filters, and reflect the un¬ 
certainty otherwise, via smaller |C,;| and larger Di 
values. 


model comparisons which are not relevant to a particular 
classification type, contribute significantly to the noise in 
the clustering process. We also attempted a clustering on 
the differences between the model LLOOCVs and AICcs, 
instead of on the signs of their differences, in a single 
clustering step, as well as in a hierarchical supervised 
method as discussed in this paper. However, we found 
that in both cases, the number of misclassifications is 
larger due to the associated variance in the values of the 
differences in LLOOCV and AICc, which is mitigated by 
reducing them to binary statistics using the sign of their 
differences alone. 

We combine the clustering classifications from the sec¬ 
ond clustering stage in each filter, by defining two mea¬ 
sures; a quality factor Ci which is the average of classifi¬ 
cations across the filters: 


C,= 


E f Cu 


TV 


filters 


( 10 ) 


and, the average clustering square distance Di across 
the filters: 


A = 


Ef Dij 

Nfilters 


( 11 ) 


Note, that it is favorable to assume more clustering 
centers in any filter than there are. For example, we 
could assume that a certain filter has 3 clustering cen¬ 
ters corresponding to BL, SV, or NM while it may so 
happen that one of the BL or SV centers is repeated, 
or two of the centers are relatively proximal, implying 
that the clustering really occurs only around two cen¬ 
ters. Post clustering, we filter out sources which have 
been clustered around NM centers in at least 3 bands 
they are detected in, and label them NM sources. The 
remaining sources then, are classified in at least 2 bands 
they are detected in as either SV or BL. To detect their 
type more precisely than just using their comparisons to 
the no-model, we construct another relative sign vector 
BLSVij comparing the fitness statistics of the BL mod¬ 
els directly to those of the OU process for each source, 
band-wise: 


BLSVij = { 

Sgn(LLOOCVGaussian,iJ - LLOOCVou,i,f) > 
sgn{LLOOCV Gamma , u - LLOOCV ou ,ij), 
Sgn{LLOOCVAnalytic,i,f — LLOOCV()U,i,f )) 
Sgn(AICCGaussian,i,f AICcou,i,f)) 
Sgn(AICcGamma,i,f AICcQjjfJ ), 
Sgn(AICCAnalytic,i,f AICcqu^j 

} (9) 

We aggregate BLSVij band-wise and perform a two- 
center K-means clustering ( I\ = 2) to segregate the BL 
and SV sources. We find that this type of hierarchical 
supervised clustering, i.e. filtering out the NM sources 
in the first stage and classifying the remaining sources 
as BL or SV in the second stage, is more efficient as 
compared to using all the model statistical comparisons 
concurrently in a single clustering step. This is because, 


BL sources have Ci closer to 1 while stochastically vari¬ 
able sources have Ci close to —1. Di is a measure of the 
overall reliability of the classification that decreases with 
increasing Dj. Therefore, sources which are purely BL 
will have C t = 1, Di = 0, while purely stochastic vari¬ 
ables will have Ci = —l,Di = 0. Intermediate values 
of Ci indicate disagreements between some of the band- 
wise classifications, while larger values of Di indicate a 
disagreement between the models in a given band. 


5.1. Tests On a Verification Set 

To test our classification method we constructed a reli¬ 
able verification set with a diverse range of SNe (BL) and 
AGN (SV) in order to capture, as much as possible, the 
full range of their time-variability properties. For AGN, 
we created a verification set from two sources: 1) 58 
UV-variability selected AGN with associations with PS1 
alerts within 1" from the GALEX Time Domain Survey 
(TDS) ( Gezari|2013 ) with no available spectroscopy, and 
2) 125 spectroscopically confirmed AG N PS1 alerts asso- 
ciated with galaxy hosts from SDSSfShen et al. 2011) 
and from a multipurpose Harvard/CfA program with 
the MMT to observe PS1 transients (PI Berger). The 
GALEX AGN were selected from UV variability at the 
5er level in at least one epoch, and then classified us¬ 
ing a combination of optical host colors and morphology, 
UV light curve characteristics, and matches to archival 
X-ray, and spectroscopic catalogs. The SN verification 
set consists of 131 spectroscopically confirmed Type-la, 
Type-Ib/c, Type-II, Type-IIn, and Type-IIP SNe from a 
combination of PS1 spectroscopic follow-up programs us¬ 
ing Ge mini, Magellan, and MMT described in [Rest et al.| 
(2013) and Berger et al. (in prep.). In order to test the 
performance and efficiency of our algorithm in the classi¬ 
fication of AGN and SNe, we have constructed a diverse 
and robust verification set that should be representative 
of these populations in our sample. 
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Fig. 8. — Densities of verification set SNe (left) and AGN (right) on the Ci vs D t plane. Since AGN classifications for D, > 4 occupy 
both the BL (Ci > 0) and SV(Ci < 0) regions, we only rely on classifications with D, <= 4. As a result, the SNe are classified with 93.89% 
completeness and 90.97% purity while the AGN are classified with 57.92% completeness and 95.00% purity. 


Fig. [8] shows the contours of Ci vs Di for the spectro¬ 
scopic SGN and SNe. The SNe cluster around the region 
Ci > 0.5 and Di < 4, while AGN predominantly occupy 
the regions defined by Ci < 0 and Di < 8. In general, 
the degree to which an object is BL as opposed to SV 
increases with Ci. Some AGN light curves may show 
bursting-type behavior resulting in their being classified 
in more than 1 filter as BL, consequently having Ci > —1. 
Also, AGN clustering are less reliably classified as evi¬ 
denced by systematically larger Di as compared to the 
SNe. 

This is possibly because, the damped random walk 
model may be over-simplifying the description of the un- 
derlyi ng AGN variabilit y. A more complex m odel such 
as in Kelly et al. (2010) or Kelly et al. (2014|), may be 
required to capture the true variability. Another source 
of confusion is that some of the AGN may resemble BL 
light curves if they fluctuate above the detection thresh¬ 
old for only a brief period during the observing baseline. 
Consequently, 47.88% of the verification set AGN clas¬ 
sified with Di > 4, indicating unreliable classifications. 
In order to maximize the purity of our classifications, 
with a sacrifice to completeness, we use Di = 4 as the 
bound for the classifications, below which 57.92% AGN 
and 93.89% SNe, are recovered with 95.00% and 90.97% 
purities respectively. It is possible to include other pho¬ 
tometric properties like color or host-galaxy properties 
to improve the completeness of the AGN classifications, 
however, since the focus of our present work is to only 
use time-variability as a tool for classification, we reserve 
this for future work. 

In multi-epoch surveys such as Pan-STARRSl, it may 
be possible to differentiate between AGN and SNe sim¬ 
ply by comparing variability between observing seasons. 
For example, an SN will have only one season in which 
the reduced x 2 (X 2 ) is larger than 1, while an AGN light 
curve can show variability in both seasons with xt >> L 
In Fig{l2]we test this simplified method by plotting the 


minimum of the seasonal xt s for the verification set AGN 
and SNe. If xt = 5.75 is used to separate AGN from SNe, 
then 55.73% of the AGN can be recovered with 86.45% 
purity. We find that, using this value of xl t° sepa¬ 
rate AGN and SNe, results in maximal purity for the 
AGN sample. While this method yields a completeness 
for AGN that is similar to that of our light-curve clas¬ 
sification algorithm, for SNe, the performance is much 
worse, with 87.69% of SNe recovered at 47.22% purity. 
This high contamination rate for SNe is due to the ex¬ 
tensive overlap between the AGN and SN in the region 
xt < 5.75. Hence, we conclude that for maximum purity, 
a more sophisticated method such as the one adopted in 
this paper is necessary. 

5.2. Final Classifications and Properties of 
Extragalactic Sources 

We begin classifying our 4361 extragalactic transient 
alerts by first selecting out sources which are clustered 
around the NM center in at least 3 of the 4 bands. We 
find 570 such sources (NM sources hereafter). Visual in¬ 
spection of the NM source light curves reveals that the 
majority are the result of noisy image-differencing light 
curves, most often due to large excursions in flux from 
image differencing artifacts, and not statistical errors due 
to faint fluxes, in most of the bands. This can be seen in 
Fig- E which shows the minimum source magnitude in 
the gpi,rpi,ipi , and zp± bands for all sources, includ¬ 
ing that for NM sources (dark green), which barring the 
brightest end of the magnitude distribution, follows the 
overall magnitude distribution of extragalactic sources 
(black), indicating no strong biases toward fainter mag¬ 
nitudes. 

Fig. [13] shows Ci vs Di contours for the 3791 extra¬ 
galactic sources classified SV and BL. We determine that 
there are 2262 SV sources and 1529 BL sources in the 
dataset. Fig. [14] shows the distributions of SV, BL, and 
NM by MD field, with SV being the most common class 
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Fig. 9.— Distribution of SV (blue), verification set AGN (cyan), BL (red), verification set SN (orange), and NM (dark green) as a 
function of minimum magnitudes in the <?pi, rpi, ipi, and zpi bands. 




TABLE 2 

Source variability and offset classifications 


Type 


Nuclear (offsets < 0.55") Off-Nuclear (offsets > 0.55") 


Burst-Like 
Stochastic Variable 
No-Model 


of extragalactic alert in all fields. 

We combine the light-curve classifications, with host 
galaxy offsets, in order to define a robust photometri¬ 
cally selected sample of AGN and SNe, from nuclear 
SV and off-nuclear BL, respectively. In order to de¬ 
termine our cut-off for off-nuclear sources, we first fit 
the offset distribution for all the sources, with a bi- 
modal distribution (Fig. 10), to derive a nuclear (predom¬ 
inantly AGN) population with n nuc = 0.26, <J nuc = 0.14, 
and off-nuclear (predominantly SNe) population with 
[J'off—nuc — 0.48, (Joff—nuc — 0.3T. SNe can be coinci¬ 
dent with galaxy nuclei due to the limited spatial res- 


689 812 (SNe) 

1233 (AGN) 1027 

449 121 

olution of the images. AGN, however, should not have 
significant offsets from their host galaxy centers, unless of 
course, they are more exotic objects such as recoiling su- 
permassive black holes, or dual AGN. Therefore, objects 
with offsets greater than /i nuc + 2 a nuc = 0.54 are most 
likely to be SNe. For offsets < 0.54” the AGN popula¬ 
tion is more than 97% complete, and there is negligible 
contamination in the supernova population by AGN, be¬ 
yond this offset. Following this line of reasoning, we use 
the AGN from our verification set to determine the nu¬ 
clear offset distribution, shown in Fig. 


11 This is fitted 


with a Gaussian and leads to a p. agn + ‘Z&agn = 0.55 
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Fig. 10. - Bimodal distribution of PS1 extragalactic alert host 
offsets, separated into nuclear, and off-nuclear distributions, with 
A»nuc, crime = 0.26, 014 and fCoff-nuc) cr 0 ff_ nuc = 0.48, 0.37. Sources 
offset from their host galaxies by more than fi n uc + 2cr n uc = 0.54” 
are predominantly SNe. 



Fig. 11.— Offset distribution of verification set AGN which is 
well approximated by a Gaussian with hagn■> &AGN = 0.25,0.15. 
This is approximat ely t he same as the distribution for nuclear off¬ 
sets obtained in Fig |10| from the bi-modal assumption for the entire 
extragalactic alerts population. 

cut-off for AGN, which we adopt for demarcating nuclear 
sources from off-nuclear ones. The offset distribution for 
each variability class is shown in Fig. [15} The distribu¬ 
tion of SV offsets is broader than that of the verification 
set AGN, however, the broader distribution likely reflects 
the larger errors in the image difference and host galaxy 
centroids for fainter AGN, that are not represented in 
the verification set. The BL distribution is seen to ex¬ 
tend well beyond the nuclear AGN distribution, as would 
be expected for SNe. There are also a significant num¬ 
ber of SVs with offsets > 0.55”, however, we find that 
these sources are faint with a mean g ~ 21, resulting in 
their offsets being poorly determined. Table [2] shows the 


Fig. 12.— Minimum seasonal xt of the spectroscopic verification 
set of AGN and SNe in the g band. A cut-off of xt > 5.75 can be 
used to demarcate AGN from SNe, albeit with a high contamina¬ 
tion rate for SNe. 

number of sources in each variability class divided into 
nuclear (offset < /r nuc + 2a nuc = 0.55") and off-nuclear 
(offset > 0.55"). We designate the nuclear SV to be 
AGN, of which there are 1233, and the off-nuclear BL 
to be SNe, of which there are 812. In the next section, 
we use these variability/offset selected AGN and SNe, to 
define photometric priors for their easy identification in 
future surveys. 


6. CHARACTERIZATION OF VARIABILITY SELECTED 
EXTRAGALACTIC AGN AND SNE 

For upcoming multi-band, multi-epoch surveys such as 
LSST, we have shown that light-curve characterization 
combined with host galaxy offsets is a robust way to se¬ 
lect AGN, SNe, and other exotic events, and does not 
require data external to the survey such as spectroscopic 
follow-up. Using all the gpi, rp±, ip±, zp\ bands offers a 
redundancy that increases the confidence of source clas¬ 
sification. With our photometrically selected samples of 
AGN and SNe, we now characterize their key observed 
source and host galaxy properties, with the hopes of find¬ 
ing priors that can accelerate their identification in future 
surveys. 

We use the *pi-band to characterize the host galaxy 
magnitudes of our sources, since the ipi-band has the 
highest signal-to-noise ratio amongst all the PS1 bands, 
and contamination of host galaxy flux by a central AGN 
is minimized as compared to the bluer bands. Fig. |16| 
shows the distribution of host galaxy ipi for AGN and 
SN. AGN host galaxies appear significantly brighter in 
the i-band than the SN host galaxies. Preliminary red- 
shift estimates of the transient alert host galaxies indi¬ 
cate that SN h ost g alaxies have a larger mean redshift 


distribution (Heinis et al.| 
AGN host galaxies, there b) 


| in prep.p as compared to the 
thereby resulting in the observational 

bias. 

AGN themselves are much fainter in difference flux 
as compared to their host galaxy flux, and we can use 
this to further separate the AGN from the SN using 
the distribution of the differences between the minimum 
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Fig. 13.— Density maps for BL (left) and SV(right) as a function of Ci and Dj. 


(b) 



Fig. 14.— Distribution of SV,BL, and NM sources across the 10 MD fields. 
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Fig. 15.— Host galaxy offset distributions for SV, BL, and NM 
sources in arcsec. Dashed line indicates the offset above which a 
source is considered “off-nuclear”. 



Host i-band (mag) 

Fig. 16.— Distribution of SN and AGN host galaxy i-band 
magnitudes. AGN host galaxies are « 3 mag brighter than SN 
host galaxies. 


source i-band difference-magnitude and the host mag¬ 
nitude (■ imin — ihost) (Fig. 17). AGN peak variability 
amplitudes are significantly fainter (« 4 mag) relative to 
their host galaxies, as compared to that for SNe (~ 2 
mag), consequently being more difficult to detect. 

We find that by using only 4ost and (z min - Zhost), we 
can compute informative priors for the source-types from 
their relative probabilities of occurrence. Figjl8] shows 
the contours of AGN and SN in ihost and i m i n — ihost 
space. Although the AGN and SN distributions overlap 
in this space, there is a clear divide between their highest 
density regions, making it possible to separate them and 
assign relative probabilities in the overlap regions. Ap¬ 
proximating and smoothing the S Ne a nd AG N, Zh nst and 
(*min — ihost) distributions (in Fig|l6|and Fig[l7] respec¬ 
tively) by Gamma distributions, we obtain their respec- 


Fig. 17.— Distribution of the differences between the minimum 
i-band difference magnitude and the host galaxy i-band magnitude 
for all source types. AGN fluxes are typically much fainter relative 
to their host galaxies with typical (£agn — ihost ) ~ 4 mag, while 
SNe are typically 3 mag fainter than their host galaxies in the 
i-band. 

tive joint probability distributions in both parameters 
as: 


Pagn = liimin - ihost , k = 25.360 ,6 = 0.230) 

x liihost, k = 4.911 ,9 = 0.651) (12) 

PSN = T {imin ihost 5 k = 12.852, 0 = 0.400) 

x liihost , k = 11.080, 9 = 0.469) (13) 

If Nagn and Nsn are the observed number of AGN 
and SN, the relative AGN likelihood for any set of values 
ihost and 'imin - i h o S t is given by 


_ NagnPagn 

Pagniagn,sn) — -jG -; y;- 

A AGNPAGN + J^SNPSN 


(14) 


where Pagn|agn,sn) is the probability that an object 
is an AGN, given that it is either an AGN or an SN. 
Assuming that the number of AGN and SNe scale lin¬ 
early with the number of SV and BL sources respec¬ 
tively, we obtain Nagn = 2262 and Nsn = 1529. 
Fig{l9] is a smoothed version of Fig{l8] and shows the 
contours of Pagniagn.sn- SNe at their brightest, being 
brighter than AGN, and also being the source-type that 
dominate alerts from fainter galaxies, typically occupy 
smaller z mira — ihost and larger ihost (redder contours). 
Whereas, AGN fluxes being smaller compared to their 
host galaxy fluxes, and from less distant galaxies, occupy 
larger imin — ihost and smaller ihost (bluer contours). The 
probability of a SN at any given point in this parame¬ 
ter space is 1 — pagn|agn,sn- The verification set AGN 
(black stars) and SNe (magenta circles) are plotted for 
reference. 

For the problem of classification in real-time from a 
large data stream such as the LSST transient alerts, we 
have found that for a magnitude-limited survey, simply 
using the i-band peak source magnitude and i-band host 
magnitude as priors, one can produce a robust prelimi¬ 
nary AGN vs. SN classification, in order to help filter out 
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a sample for more tedious methods such as spectroscopic 
or time-series identification of sources. 


0 20 40 60 80 

AGN - % of maximum density (color line contours) 



Host i (mag) 


0 20 40 60 80 

SN - % of maximum density (gray shaded contours) 

Fig. 18. — Distribution of AGN (contours in blue to red) and 
SN (contours in gray to black) ihost and i m i n — ihoat • A clear 
separation can be seen between the highest density regions of the 
two source types. 


7. CONCLUSIONS 

In this paper, we discussed a multi-band difference- 
flux time-series based method for the classification of 
4361 PS1 MD extragalactic difference-image sources into 
stochastic and bursting. Using a star-galaxy catalog to 
select extragalactic sources, we classify them into SV and 
BL sources using band-wise difference-flux characteriza¬ 
tion. Although this method can use actual or difference- 
magnitude time-series, difference-flux time-series are pre¬ 
ferred over difference magnitudes which are log scaled, 
thus circumventing the problem of negative difference- 
flux excursions in SV light curves, for which magnitudes 
cannot be defined. We use multiple BL models to model 
the shapes of BL light curves, an OU process to model 
SV light curves, and a No-Model to identify white-noise 
dominated light curves. Since the models only attempt 
to differentiate between coherent single-burst type be¬ 
havior and stochastic variability, they do not assume any 
underlying physical processes for the sources, making the 
method widely applicable. The use of multiple BL mod¬ 
els is justified for statistical redundancy in the param- 
eterizations of the light curves, as well as for modeling 
the gamut of shapes of BL light curves. We find that 
using all 3 BL models in the place of any one or two 
of them, improves the purity and completeness of the 
variability classifications of our AGN and SN verifica¬ 
tion sets. We estimate the model fitnesses using their es¬ 
timated corrected-Akaike information criteria, and their 
leave-out-one-cross-validation likelihoods in each filter. 
The use of these independent derived statistical mea¬ 
sures, one of which is suited to simply assess light curve 
shape characteristics, and the other to assess the over¬ 
all robustness of the model, works to fortify the derived 
classifications. 

We then construct decision vectors RVij for each 


source, based on the AICc and the logarithm of the 
LOOCV for all the time-series models, which are com¬ 
bined in two clustering steps across the sources, and clas¬ 
sified using a supervised K-means clustering method to 
arrive at the final filter-wise classifications; we filter out 
the NM sources in the first step, and then we separate out 
the SV and BL sources in the second. The use of time- 
series in multiple bands increases the reliability of our 
classifications. We then define two quality measures C) 
and D t , which are filter-wise averages of the final clus¬ 
tering classification parameters, in which space the SV 
and BL can be separated. We find that our method re¬ 
sults in 183 verification set AGN being classified with 
95.00% purity and 57.92% completeness, and 130 verifi¬ 
cation set SNe classified with 90.97% purity and 93.89% 
completeness. We use our method to classify all the ex¬ 
tragalactic difference-detection alerts into 2262 SV, 1529 
BL, and 570 NM best-fit sources. We then construct a 
robust photometrically selected sample of 812 SNe and 
1233 AGN, using a combination of light-curve class and 
host galaxy offset, to characterize their variability prop¬ 
erties and the properties of their host galaxies, in order 
to better inform their real-time identification in photo¬ 
metric surveys. We find that simply the combination of 
i-band host magnitude and the i-band source magnitude 
can be used as robust preliminary indicators of source 
type. 

We demonstrate that our method can be used to sep¬ 
arate SV from BL using the self-contained data (multi¬ 
epoch image-differencing and deep stacks) available in 
multi-band time domain surveys, such as PS1 and LSST. 
However, one could go further and use other variabil¬ 
ity based parameters in conjunction with our time- 
series method, together with host galaxy offsets, col¬ 
ors, and morphology, and external information from 
multi-wavelength catalog associations, in a larger, com¬ 
prehensive hierarchical classification scheme to improve 
classification accuracy, characterize known sub-classes of 
sources, as well discover new classes of sources. In addi¬ 
tion to the classification of variables and transients into 
broad general classes and particular sub-classes via the 
use of exact models, ensemble studies of their general 
properties can be readily performed; for example, the 
general properties of the host galaxies of AGN and SNe; 
the rates and properties of SNe and their subclasses; the 
variability timescales and amplitudes of AGN and their 
subclasses, and subsequently, the estimation of the black 
hole mass function; are some of the questions that can be 
readily answered using the model-fit parameter distribu¬ 
tions for the respective classes. In the era of wide-held 
synoptic surveys generating millions of transient alerts 
per night, such self-contained photometric identification, 
classification, and characterization of transients based on 
light-curve characteristics and host galaxy properties will 
be essential. 
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'host < ma 3) 


Fig. 19. — Smoothed distribution of relative AGN probability (Eq |14[ ) in the i^ost - imin ~ *host plane. The probability distributions, 
derived from the density of the photometrically selected AGN and SIN samples in each parameter, are smoothed and approximated by 
Gamma distributions. The overall distribution is obtained by multiplying the distributions in each parameter. The verification-set AGN 
(black stars) and SNe (magenta circles) are plotted for reference. 
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